Finite- Temperature Monte Carlo Calculations For Systems With Fermions 



Or 



On 
O 
On 
O 
OS 

o> 

-i— > 
c3 



O 

o 



X 



Shiwei Zhang 

Department of Applied Science and Department of Physics, College of William and Mary, Williamsburg, VA 23187 

(February 1, 2008) 

We present a quantum Monte Carlo method which allows calculations on many-fermion systems 
at finite temperatures without any sign decay. This enables simulations of the grand-canonical 
ensemble at large system sizes and low temperatures. Both diagonal and off-diagonal expectations 
can be computed straightforwardly. The sign decay is eliminated by a constraint on the fermion 
determinant. The algorithm is approximate. Tests on the Hubbard model show that accurate results 
on the energy and correlation functions can be obtained. 
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The quantum Monte Carlo method for simulating 
grand-canonical ensembles, originally formulated by 
Blankenbecler, Scalapino, and Sugar (BSS) 0], is widely 
applied in areas spanning condensed-matter, high-energy, 
and nuclear physics. The method allows essentially exact 
calculations of finite-temperature equilibrium properties 
of interacting fermion systems. It expresses the partition 
function as a many-dimension integral over a set of ran- 
dom auxiliary fields. The many-dimensional integral is 
then computed by Monte Carlo (MC) techniques. 

As all current fermion quantum Monte Carlo methods, 
however, the BSS algorithm suffers from the well-known 
sign problem |^|,^| . The integrand of the partition func- 
tion is not all positive. Indeed its average sign approaches 
zero as the temperature is lowered. As a result, contribu- 
tions from the Monte Carlo samples largely cancel. The 
partition function, which is given by the difference, be- 
comes a vanishingly small quantity compared to the MC 
noise. The computational cost for fixed statistical ac- 
curacy scales exponentially with system size and inverse 
temperature. While for many problems the BSS algo- 
rithm is the most, sometimes only, feasible approach, 
the sign problem has remained completely uncontrolled 
in the algorithm. This has severely limited the temper- 
atures and sizes accessible, and has prohibited studies of 
a variety of interesting problems in correlated systems, 
particularly concerning true phase transitions. 

In this Letter, we present a finite-temperature method 
which is free of any decay of the average sign and which 
retains many of the advantages of the BSS formalism, 
thus allowing grand-canonical calculations at lower tem- 
peratures and larger system sizes with favorable scaling. 
Below we first derive a set of exact constraints on the 
auxiliary fields which eliminates any negative contribu- 
tion to the partition function. An approximation is then 
made to impose these constraints in the MC sampling to 
control the sign problem. We develop an algorithm to 
effectively carry out the MC sampling under the approx- 
imate formalism. We illustrate the method by applying 
it to the one-band Hubbard model. We show that accu- 
rate results, on both the energy and various correlation 



functions, can be obtained with the new method, even 
with simple forms of the approximate constraint. 
The expectation value of a physical observable O is: 
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where ft = 1/ kT is the inverse temperature. The chemi- 
cal potential term is implicit in the Hamiltonian H. The 
partition function in the denominator can be written as 



Z = Tr(e-P H ) = Tr[e 



-AtH 



, e -ArH e -Ar H] 



(2) 



where At = /3/L and L is the number of "time slices" 
on the right-hand side. 



We next write the many-body operator e~ ArH in terms 
of single-particle operators. This is possible for most 
Hamiltonians or Euclidean actions of interest. For ex- 
ample, the Hubbard-Stratanovic transformation 0] can 
be applied for a Hamiltonian H which contains one- 
and two-body terms, denoted by K and V, respectively. 
This transformation replaces the two-body term e~ ArV 
by one-body interactions with a set of random exter- 
nal fields. Combining the result with the one-body term 
e -ArK we can wr jt e 



-AtH 



(3) 



where x denotes the random external auxiliary fields and 
-B(x) is a single-particle operator. The sum over all aux- 
iliary fields recovers the interaction. For simplicity we 
have written the integration over x as a discrete sum. 
We have also suppressed spin indices, as well as the dis- 
tribution function of x. The approximation in Eq. (^) is 
from the Trotter error, which is of C(Ar 2 ) or less. 

In the standard BSS formalism, Eq. (||) is substituted 
into Eq. (||). The trace over fermion degrees of freedom 
is then performed analytically IllJq], which yields 



Tr 



= ]Tdct[/ + B(x L ) • • • B(x 2 )B(xi)], 
x 
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where X = {xi,X2, ■ ■■ denotes a complete path in 
auxiliary-field space. If the size of the single-particle ba- 
sis (e.g., number of spatial lattice sites) is N, the single- 
particle propagator -B(xj) is an N x N matrix and / is 
the corresponding unit matrix. The fcrmion determinant, 
which we will denote by D(X), can be computed for each 
X. The sum over all paths can therefore be evaluated by 
MC methods. However, D(X) is not always positive. 
As illustrated in Fig. [l]a, the physical contribution comes 
from the small difference between the positive and nega- 
tive components. The MC samples of X are drawn from 
the probability distribution defined by |Z?(A)|. As (3 in- 
creases, D(X) approaches an antisymmetric function and 
its average sign vanishes exponentially. The variance in 
the MC estimate of Eq. (|l|) thus diverges, and the sign 
problem occurs. 

A main obstacle to understanding and controlling the 
problem lies in the implicit and complex nature of the 
path-integral picture in this formalism. To gain insight, 
we return to the original form of Z in Eq. (|^). We 
will use B to denote e~ ArH and imagine the following 
thought experiment to generate all possible auxiliary- 
field paths X. Beginning with Tr[BB ■ ■ ■ BB], we sub- 
stitute B with Eq. (||), one at a time from right to left. 
After I such steps, the partition function can be written 
as E{x 1 ,x 2 ,...,x ! } 7 'K{ x i. x 2,- ' ' ,3Q},B), where P t is 

P i ({x 1 ,x 2 ,---,x i },6) 

= Ti[BB-_BB( Xl )---B( X2 )B(x 1 )}. (5) 

L-l 

As we proceed, we construct paths by including all pos- 
sible values of x;. After L steps, all iTs are replaced and 
all complete paths X are generated. Note that, while not 
the case in general, the trace in Eq. (||) can be performed 
when I — L, which, as expected, gives D(X) of Eq. (||). 

We now examine the procedure more closely, first at 
At — > 0, where Vi is continuous in I, the length of the 
partial path. In particular, we consider the case when Vi 
becomes zero for a certain partial path {xi,X2, • • • , x;}. 
This means that, after the remaining L — I steps have 
been finished, the sum over all possible configurations 
of {jq+i,jq+2j • • • ,xl} will simply reproduce the i3's in 
(|J), leading to zero by definition. In other words, any 
complete path whose first / elements are {xi, x 2 , • • • ,X;} 
is "noise"; the contributions of such paths cancel in Z. 
The signature of a noise path is Vi = for at least one 
I. Since Vq > 0, this shows that a complete path con- 
tributes if and only if the following L conditions hold: 

Pi({xi,X2,---,Xi},B)>0, ? = 1,2,--.,L. (6) 

If we impose the constraints in Eq. (^|) in our proce- 
dure to generate the paths, we can eliminate all noise 
paths while selecting all contributing paths. The con- 
straints are equivalent to having an absorbing boundary 



at the V\ = axis in Fig. lb, thereby making the proba- 
bility distribution of the generated complete paths van- 
ish smoothly at the axis. This boundary condition (BC) 
eliminates complete paths that come in contact with the 
axis at any point, which cancels out the antisymmetric 
part of D(X) in Fig. la. The algorithm remains exact. 




FIG. 1. Schematic illustration of the sign problem and the 
constraints to control it. Fig. (a) shows the integrand D(X) of 
the partition function Z. The X-axis represents an abstrac- 
tion of the many-dimensional auxiliary- field paths X; each 
point denotes a collection of X's, e.g., in the sense of a bin 
in a histogram. In standard MC, |D(X)| is sampled, while 
only the shaded area contributes. Fig. (b) shows Vi (Eq. (^)) 
as a function of the length of the partial path, I, for several 
paths. When Vi becomes 0, ensuing paths (dashed lines) can- 
cel. Only complete paths with Vi > for all I (solid line) 
contribute in Z; they lead to the shaded area in (a). 

At finite At, paths are discrete. But the BC is the 
same for the underlying continuous paths. To the lowest 
order, the constraints in Eq. (^|) allow imposition of the 
BC under the discrete representation; the contact point 
("triple point" in Fig. lb) is approximated by the first 
I for which Vi < 0. A higher order approach, which we 
use, is to interpolate between this I and I — X, with the 
probability to terminate at I — 1 approaching 1 smoothly 
if T^z — i — > |^] . It is important to note that, in both 
approaches, the finite-Ar error vanishes as At — > 0. 

B is not known in practice. We replace it by a known 
trial propagator Bt- The constraints now yield approxi- 
mate results, which become exact if Bt is exact. If Bt is 
in the form of a single-particle propagator, we can analyt- 
ically evaluate the trace in Eq. (||) by making use of the 
same identity that produced Eq. (Q). The constraints 
in (J^j can now be written as: 

L-l 

Vj = dot [/+([] B T ) B( Xl ) ■ ■ ■ B(X!)] >0 (7) 

TCI— 1 

for each / on 1 < ( < I, where we have introduced the 
shorthand Vj for Vi({xi, x 2 , ■ • ■ ,Xj}, B T ). 

The idea of the new method is then to generate MC 
samples of X which both satisfy the conditions in (ffl) and 
are distributed according to D(X) Q. To realize this 
efficiently, we construct the following algorithm, which 
builds directly into the sampling process both the con- 
straints and some knowledge of the projected future con- 
tribution. In terms of the partial contributions Vj , the 
fcrmion determinant D(X) can be written as 
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We construct the path X in L steps, corresponding to 
stochastic representations of the L ratios in Eq. (@). Wc 
start from Vq, i.e., L Bt's in place of i?'s, with overall 
weight 1. Then, successively from I = 1 to L, we: (a) 
pick an x; from the conditional probability density func- 
tion p(xj|x l _ 1 ,-.-,x 2 ,xi) > defined by (Vf/VfLJ/C 
and, (b) multiply the overall weight by the normalization 
factor C = I^T-v The algorithm allows x; to be 

selected according to the best estimate of its potential 
contribution, reflecting the integrated (i.e., with dashed- 
line paths in Fig. (lb) already canceled out) effect of 
all subsequent paths from xj. Note that the probability 
distribution for x; vanishes smoothly as Vj approaches 
zero, and the constraints are naturally imposed. 

We simultaneously propagate an ensemble of paths. 
The contribution of each path X in Z is given by its fi- 
nal weight. Given X, we can calculate both equal-time 
and (imaginary) time-dependent correlations through the 
single-particle Green's functions @]. The expectation in 
Eq. (|l|) is a weighted average over X. The statistical ac- 
curacy improves as the procedure is repeated and more 
paths are generated. 

We mention several technical issues, (i). We have 
chosen a non-interacting propagator, e~ ArK , as Bt- 
More general mean-field propagators, including ones 
with imaginary-time dependence, can be incorporated 
straightforwardly, (ii) . We divide each step for each path 
into sub-steps, in which we apply (a) and (b) to individ- 
ual components of x;. This simplifies p and C (of the 
sub-steps) Q,|(|. (iii). As paths are evolved, products of 
-B(x) and Bt must be stabilized against round-off errors 
H . (iv). The weights of paths fluctuate as they are prop- 
agated. We apply a population control mechanism |^] to 
improve efficiency. A detailed account of these and other 
algorithmic issues will be published elsewhere. 

The algorithm we have described provides the finite- 
temperature counterpart of the ground-state constrained 
path Monte Carlo (CPMC) method ||. The latter, which 
has been applied to study various lattice models, elimi- 
nated the sign decay in T = OK auxiliary-field calcula- 
tions by constraining paths in Slater determinant space 
with a trial ground-state wave function \tpr) The 
chief difficulty in generalizing the concept of a constrain- 
ing wave function or density matrix fTl|| to the finite- 
temperature formalism is two-fold: (i) In this formalism, 
paths do not originate or end at the same point in Slater 
determinant space; different paths would thus require dif- 
ferent constraining conditions. Indeed paths do not even 
have the same "dimension" . (ii) With the analytical eval- 
uation of the trace, the path-integral picture is implicit 
and would likely prevent implementation of such con- 
straints. The new algorithm overcame the difficulty It 
also provides a unified view of the zero- and finite- T algo- 



rithms. The constraining \iPt) in T = OK CPMC can be 
understood in terms of Bt operating on an initial state. 

We now apply the new algorithm to study the one-band 
Hubbard model. The model consists of interacting elec- 
trons on a square lattice. The Hamiltonian H = K +V is 
given by K = -t Y,{ij)a( c L c j<* + h.c.) - MEi( n iT + n il) 
and V — UJ2i n i1 n ili where c\ a creates an electron of 



spin a on site i, ni a — c\ a Ci a is the electron number 
operator, and ( ) indicates near-neighbors. The on-site 
Coulomb repulsion is U > 0. In connection with high- 
T c superconductivity, the Hubbard model has been the 
subject of intense theoretical effort for the past decade. 
The model provides a good test case, with both its chal- 
lenging nature and the availability of certain benchmark 
data. Quantities of particular theoretical and experimen- 
tal interest include the momentum distribution n(k) and 
the (i-wave electron pairing correlation Pd (1) flj] . 

We study lattices of size \A/V x y/~N with periodic 
boundary conditions. The desired electron density (n) = 
(J2ia n i<r) /N is achieved by adjusting \i. Our trial prop- 
agator Bt is e~ ArK multiplied by e ~ ArVT 12 iCT n ^ f where 
vt is a parameter. The second term in Bt accounts for 
e -Ar\/ j n j.jjg sense f res tricted Hartree-Fock. 

In Fig. H and in Table I, we show results for a 4 x 4, 
U — 4 system where the sign problem is the most severe. 
This limits the range of temperatures where accurate cal- 
culations can be done with the standard algorithm. At 
[3 = 12, the average sign in BSS, (s), is projected to be 
less than 0.01 from the exponential decay rate H and the 
numbers in Table I; this /3 is thus not reachable by BSS 
with present computing power ^M. The system hence 
presents a challenging test case for the current algorithm. 
At high T, our algorithm gives results in excellent agree- 
ment with BSS results |L5|, which are exact. At low T, it 
reaches convergence and leads to results consistent with 
those from ground-state CPMC and in good agreement 
with those from T = K exact diagonalization [|l6| . 

In Fig. ||, we show new results for an 8 x 8 lattice. The 
electron filling of (n) = 0.82, which is in the physically 
relevant region, shows the worst sign problem, with (s) in 
BSS falling to ~ 0.1 at (3 = 6 ||. Accurate and system- 
atic calculations have therefore not been possible on this 
system. The new algorithm, on the other hand, required 
only modest computing time (about 2 days on a single 
processor of an SGI Origin200 workstation for (3 = 16) 
to reach the excellent statistical precision shown in the 
figure. As T decreases, the Fermi surface appears to con- 
tract along (it, 7r), while bulging along (tt, 0). The d-wave 
electron pairing correlation at large pair separations in- 
creases with decreasing T. The non-interacting system, 
however, also shows the same behavior. In fact, -Pd(|l|) 
in the latter is larger than the corresponding interacting 
results, consistent with observations from ground-state 
CPMC |17|. More systematic calculations, at different 
(n), U, and system size, are currently being performed. 
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In summary, we have presented a quantum MC algo- 
rithm which allows finite-temperature, grand-canonical- 
ensemble simulations of fermion systems without any 
decay of sign. The method is approximate. We have 
shown that accurate results can be obtained with a sim- 
ple constraining propagator Bt- An improved Bt will 
lead to improved results, and the method becomes ex- 
act when Bt is exact. The algorithm makes possible 
calculations under the field-theoretical formalism whose 
required computer time scales algebraically, rather than 
exponentially, with inverse temperature and system size. 
With the second-quantized representation, it comple- 
ments the restricted path-integral MC method |pj| . The 
algorithm automatically accounts for particle permuta- 
tions and allows easy computations of both diagonal 
and off-diagonal expectations, as well as imaginary-time 
correlations. We expect the method and the concept 
brought forth here to see many applications, and to sig- 
nificantly enhance the applicability of quantum simula- 
tions in interacting lattice fermion systems. 
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R. L. Sugar, and J. W. Wilkins for helpful conversations, 
and the INT at the U. of Washington for hospitality, 
where part of the work was done. This work was sup- 
ported by the NSF under grant DMR-9734041 and by an 
award from Research Corporation. 
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FIG. 2. Comparison with available data for a 4 x 4 system 
with U = 4 and (n) = 0.875. The main graph shows the 
energy. The diamond at T — is from exact diagonalization. 
The inset shows the density-density correlation function be- 
tween near-neighbor sites. The algorithm accurately predicts 
the development of strong antiferromagnetic correlation as T 
decreases, despite the use of a constraining propagator Bt 
which by itself gives incorrect physics (flat line). At low T, 
the results converge to that of T = OK CPMC (triangle). Er- 
ror bars in "current" are smaller than symbol size, and are 
not shown jl3| . BSS results are from Ref. JljJ. For compar- 
ison, squares at T — 0.1667 show BSS results with the sign 
neglected, which is an uncontrolled approximation |?J. 
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FIG. 3. Temperature dependence of the momentum distri- 
bution (main graph) and d-wave pairing correlation (inset) for 
an 8 x 8 system with U = 4 and (n) = 0.82. (Recall n(k) is 
the Fourier transform of G(l).) As temperature (1//3) lowers, 
the momentum distribution, shown along two directions in 
k-space, becomes more anisotropic, and the long-range part 
of the d-wave pairing correlation increases. 



TABLE I. Further comparison of the current method with 
BSS and exact diagonalization (ED), on the same system 
as that of Fig. 2. G(l) is the average Green's function 
(ct, lCT Ci CT ), and Pd(l) the d-wave pairing correlation, at sepa- 
ration 1 = (l x , ly). The average sign in BSS is given by (s). In 



the last row, ED results are shown for G, while ground-state 
CPMC result for P^, the latter is not exact. Numbers in 
parentheses indicate statistical errors in the last digit. 
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G(1,0) 


G(2,2) 


ft(2,l) 


3 


current 




0.1631(1) 


-0.0415(1) 


0.0625(2) 




BSS 


0.99 


0.1631(1) 


-0.0418(1) 


0.0630(3) 
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current 




0.1663(3) 


-0.0470(4) 


0.077(2) 




BSS 


0.44 


0.1662(2) 


-0.0465(2) 


0.083(3) 
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current 




0.166(1) 


-0.050(1) 


0.078(2) 
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exact 




0.167 


-0.051 


0.078(2) 
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